STABILITY  RESULTS  OF  A 
CONTINUOUS  STIRRED  TANK  REACTOR 


By 

SUBRAMANIAM  PUSHPAVANAM 


A DISSERTATION  PRESENTED  TO  THE  GRADUATE  SCHOOL 
OF  THE  UNIVERSITY  OF  FLORIDA  IN  PARTIAL  FULFILLMENT 
OF  THE  REQUIREMENTS  FOR  THE  DEGREE  OF 
DOCTOR  OF  PHILOSOPHY 

UNIVERSITY  OF  FLORIDA 


1989 


ACKNOWLEDGMENTS 


The  successful  completion  of  this  dissertation  would  have  been 
impossible  if  not  for  the  assistance  and  help  I received  from  numerous 
people.  First  mention  must  be  made  of  the  guidance  and  sound  technical 
advice  I received  from  Dr.  R.  Narayanan,  my  graduate  advisor.  He  gave 
me  complete  independence  and  freedom  to  work  on  problems  of  my  choice. 
This  made  working  under  him  a real  pleasure.  The  one  other  person  who 
made  a significant  contribution  to  the  dissertation  was  Dr.  Lewis 
Johns.  I will  fondly  remember  the  motivating  and  inspiring  discussions 
I have  had  with  him.  They  would  perk  me  up  whenever  I felt  disillu- 
sioned. He  helped  focus  my  attention  on  my  dissertation  and  its  comple- 
tion. I am  thankful  to  the  other  members  of  my  supervisory  committee, 
Dr.  David  W.  Mikolaitis,  Dr.  Ulrich  H.  Kurzweg  and  Dr.  Gerasimos  K. 
Lyberatos,  for  their  time  and  effort,  and  to  the  Chemical  Engineering 
and  Mathematics  Departments  at  the  University  of  Florida  for  financial 
assistance. 

I would  also  like  to  thank  my  mother,  Yogambal,  father,  Subra- 
maniam,  brother,  Prakash,  and  sister,  Rekha,  for  their  strength  and 
support;  Dan  O'Neil,  Keith  Johnson,  Mark  Stephens,  Ken  Marshall  and 
Arunan  Nadarajah  in  the  department  for  their  friendship  and  making  me 
feel  at  home  in  this  "foreign"  country;  Raju,  my  roommate  and  friend  of 
nine  years,  for  general  moral  support  and  especially  for  putting  up  with 
my  occasional  bad  moods;  office  mate  Bill  Harter  for  his  patience  in 


11 


making  me  "computer  literate"  and  his  assistance  in  plotting  all  the 
figures  in  this  dissertation;  Preraa  of  Jacksonville  at  whose  house  I 
have  spent  many  a weekend  enabling  me  to  recover  from  the  stress  of 
graduate  school.  Two  secretaries  in  the  department  deserve  special 
mention.  Not  only  did  Debbie  Hitt  type  this  dissertation,  but  she  along 
with  Shirley  Kelly  made  life  pleasant  in  the  department.  I finally  wish 
to  thank  myself  for  doing  this  work. 


iii 


TABLE  OF  CONTENTS 


page 

ACKNOWLEDGMENTS  ii 

NOTATION  V 

ABSTRACT  viii 

CHAPTERS 

1 INTRODUCTION  1 

2 MODEL  AND  STEADY-STATE  BEHAVIOR  5 

3 LINEAR  STABILITY  ANALYSIS  11 

4 THE  D-PARTITION  METHOD  AND  APPLICATIONS  TO  THE 

CSTR  MODEL  22 

5 CONCLUSIONS  AND  SCOPE  36 

APPENDICES 

A DERIVATION  OF  THE  GOVERNING  LUMPED  EQUATIONS  OF  A 

CSTR  40 

B VERIFICATION  OF  THE  MODEL  WITH  EXPERIMENTS  IN  THE 

LITERATURE  43 

C EQUIVALENCE  OF  THE  EQUATIONS  GOVERNING  A CSTR  AND 

THOSE  OF  A SURFACE  REACTION  ON  A WIRE  45 

REFERENCES  48 

BIOGRAPHICAL  SKETCH  50 


IV 


NOTATION 


heat  transfer  area  between  reactor  and  coolant 
dimensionless  heat  of  reaction  parameter 
conversion  in  reactor 
reactant  concentration  at  inlet 

specific  heats  of  reactant  and  coolant  respectively 

steady-state  conversion 

determinant  of  matrix  M 

Damkohler  number,  dimensionless  parameter 

activation  energy  of  reaction 

heat  of  reaction 

reaction  rate  constant  evaluated  at  T^^ 
thermal  conductivity 

length  of  cooling  coil  or  catalytic  wire 
dimensionless  number,  similar  to  Lewis  number,  defined 
equation  2.3 

Jacobian  matrix  in  equation  3.1 

degree  of  characteristic  equation 

sum  of  principal  2x2  minors  of  M 

flow  rates  of  reactant  and  coolant  respectively 

universal  gas  constant 

% 

characteristic  equation  variable,  eigenvalue 
perimeter  of  cooling  coil 


dimensionless  reactant  temperature 
dimensionless  coolant  temperature 
dimensionless  inlet  coolant  temperature 
inlet  reactant  temperature 

steady-state  values  of  reactant  and  coolant  temperatures 

overall  heat  transfer  coefficient 

volume  of  reactor,  cooling  jacket  respectively 

velocity  of  coolant 

imaginary  part  of  s 

variable  defined  in  equation  3.1 

steady-state  conversion,  equals 

Greek 


dimensionless  parameter 
modified  heat  transfer  coefficient 

dimensionless  heat  transfer  coefficient,  identical  with 
Poore's  definition 

discriminant  of  the  cubic  characteristic  equation 
parameter  in  characteristic  equation 
parameter  in  characteristic  equation 
densities  of  reactant,  coolant  respectively 
real  part  of  eigenvalue  s 
determinant  defined  in  equation  ^.6 


Subscripts 


coolant  variable 
variable  at  inlet 

gas  phase  temperature  for  the  reaction  on  a wire 

Superscripts 

deviation  variable 
amplitude  of  deviation  variable 
dimensional  variables 
steady-state  variable 


catalytic  wire  variable 
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Instabilities  of  many  different  kinds  occur  in  chemical  reactor 
systems.  These  are  primarily  due  to  the  nonlinear  effects  in  the  pro- 
cesses occurring  in  these  systems.  The  simplest  of  the  systems  which 
give  rise  to  the  interesting  phenomena  as  limit  cycles,  multiple  steady 
states  are  the  first  order  reaction  in  a continuous  stirred  tank  reactor 
(CSTR).  The  strong  nonlinear  dependence  of  the  reaction  rate  on  temper- 
ature gives  rise  to  these  phenomena.  In  order  to  achieve  tractability 
of  the  problem,  it  is  common  practice  to  make  assumptions.  These  sim- 
plify the  equations  and  yet  retain  the  essential  features  of  the 
system. 

This  dissertation  can  be  divided  into  two  parts.  In  the  first 
part,  we  modify  the  equations  governing  a CSTR  by  considering  the  finite 
thermal  capacity  of  the  cooling  jacket.  This  relaxes  the  assumption 
that  the  coolant  temperature  is  a constant  at  its  inlet  value.  The 
energy  balance  equation  for  the  coolant  is  considered  and  we  study  its 
influence  on  the  stability  of  the  steady-states. 

viii 


In  the  latter  part,  a completely  new  and  different  approach  to  the 
CSTR  problem  is  discussed.  The  standard  method  of  analysis  for  the  CSTR 
is  to  fix  the  parameters  and  then  evaluate  the  steady-state  conversion 
and  its  stability.  The  limitation  of  this  method  is  we  do  not  know  the 
conversion  and  its  stability  until  we  finish  solving  the  problem.  The 
new  method  involves  calculating  the  entire  combination  of  parameter 
values  which  yield  an  asymptotically  stable  conversion  level.  The  D- 
partition  method  is  used  to  accomplish  this. 

In  the  conclusions  section,  the  evolution  of  the  dissertation  in 
its  final  form  is  discussed.  Promising  avenues  for  future  work  are 
described,  i.e.,  studying  chaotic  dynamics,  quasiperiodic  oscilla- 
tions. These  should  be  feasible  extensions  of  the  model  we  have  stud- 
ied . 
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CHAPTER  1 
INTRODUCTION 


The  exothermic  first  order  reaction 

A -»•  Y 

in  a continuous  stirred  tank  reactor  (CSTR)  has  been  studied  extensively 
in  the  chemical  engineering  literature  (Poore,  1973).  In  a typical 
reactor  (Figure  1.1)  the  reactant  A enters  the  reactor  at  a fixed  con- 
centration and  temperature.  Under  the  conditions  of  perfect  mixing  the 
effluent  concentration  and  temperature  are  the  same  as  those  existing 
inside  the  reactor.  In  a non-adiabatic  reactor  a coolant  fluid  is 
usually  circulated  through  a cooling  coil/jacket  encasing  the  reactor. 
The  coolant  absorbs  a part  of  the  heat  released  by  the  exothermic  reac- 
tion preventing  the  temperature  from  increasing  in  an  uncontrolled 
manner . 

The  strong  nonlinear  dependence  of  the  reaction  rate  on  temperature 
gives  rise  to  many  interesting  dynamic  behaviors.  The  first  elaborate 
treatment  of  this  problem  was  given  by  Poore  (1973).  He  studied  the 
dynamics  of  the  system,  obtaining  conditions  for  occurrence  of  multiple 
steady  states,  hysteresis  loops  and  limit  cycles.  His  model  contains 
three  dimensionless  parameters.  He  classified  his  3-parameter  space 
into  different  regions,  each  region  having  its  own  characteristic  behav- 
ior of  solutions.  The  algebraic  tractability  of  this  problem  has  made 
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Figure  1,1.  A typical  continuous  stirred  tank  reactor 
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it  a worthwhile  effort  for  many  investigators  to  study  this  problem. 
Most  studies  since  then  have  focused  on  (a)  studying  the  effects  of 
different  kinetics,  typically  considering  Arrhenius  kinetics  with  a 
finite  activation  energy;  (b)  using  a different  parameter,  i.e.,  flow 
rate  as  a control  variable.  The  works  of  Uppal  et  al . (1974,  1976)  and 
Balakotaiah  et  al.  (1981)  come  to  mind. 

All  analyses  of  the  non-adiabatlc  reactor  have  so  far  neglected  the 
rise  in  temperature  of  the  coolant  fluid.  This  assumption  of  a constant 
or  negligible  change  of  the  coolant  temperature  is  valid  in  the  limit  of 
infinitely  fast  flow  of  the  cooling  fluid  or  in  the  limit  of  infinite 
volume  of  the  cooling  coil/jacket.  In  the  former  case,  the  coolant  does 
not  spend  enough  time  in  the  jacket  to  gain  enough  heat  to  register  a 
significant  change  in  its  temperature.  Here  the  constant  coolant  tem- 
perature is  equal  to  its  inlet  value.  In  the  latter  case,  the  thermal 
inertia  of  the  cooling  jacket  is  very  high.  Hence,  here  once  again 
although  the  coolant  absorbs  heat  from  the  reactor,  its  temperature 
change  is  negligible  since  the  jacket  acts  as  an  infinite  sink.  Now  the 
coolant  temperature  remains  constant  at  its  initial  value. 

In  the  first  part  of  this  dissertation,  the  dynamics  of  a CSTR 
supporting  a first  order  reaction  are  modelled  by  taking  into  account 
the  effect  of  the  heat  absorbed  by  the  coolant  fluid,  i.e.,  considering 
a finite  volume  of  the  cooling  jacket/coil  and  a finite  rate  of  the  flow 
of  the  coolant.  Under  these  conditions,  the  assumption  of  constant 
coolant  temperature  breaks  down.  In  reality  the  coolant  temperature 
varies  along  the  coil.  However  as  a first  approximation  (see  Appendix 
A),  a simplistic  model  would  be  to  study  the  dynamics  of  the  coolant  by 
assuming  it  to  be  well  mixed  and  doing  an  overall  heat  balance  on  the 
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coolant  fluid.  This  enables  us  to  neglect  spatial  variations  in  the 
equation  for  the  coolant  fluid.  We  then  have  an  overall  lumped  param- 
eter model  for  our  system. 

A brief  discussion  of  the  remaining  sections  now  follows.  In  the 
next  section  of  the  dissertation,  the  mathematical  model  of  our  system 
and  the  calculation  of  its  steady-state  behavior  are  discussed.  As  a 
first  step  towards  looking  for  any  dynamical  behavior,  a linear  stabil- 
ity analysis  was  then  performed  to  obtain  conditions  for  the  onset  of 
limit  cycles.  This  gives  us  some  insight  into  the  dynamic  behavior  of 
our  model.  We  study  the  stability  and  bifurcation  of  the  steady-state 
in  terms  of  our  parameters.  The  main  objective  of  this  part  is  to 
investigate  how  inclusion  of  the  finite  thermal  capacity  of  the  coolant 
alters  the  stability  of  steady-states.  The  focus  is  on  finding  whether 
existing  unstable  steady-states  can  be  stabilized  by  changing  the  cool- 
ant residence  time. 

In  the  second  part  of  the  dissertation,  a different  approach  to 
solving  the  CSTR  problem  is  described.  The  standard  method  used  in  the 
literature  so  far  has  been  to  study  the  stability  of  the  steady  state 
conversion  for  a fixed  combination  of  parameter  values.  We  try  to 
answer  the  question  What  are  the  total  possible  combinations  of  param- 
eters which  yield  an  asymptotically  stable  steady  conversion  level? 
This  is  done  using  the  D-partition  method.  Its  applications  to  the 
standard  CSTR  (Poore,  1973)  model  and  our  model  is  discussed.  The 
usefulness  of  this  method  over  the  standard  method  is  explained.  In  the 
conclusions  section,  possible  avenues  for  future  research  are  dis- 
cussed. These  should  be  feasible  extensions  of  the  model  presented 


here. 


CHAPTER  2 

MODEL  AND  STEADY-STATE  BEHAVIOR 


The  equations  modelling  the  standard  non-adiabatic  CSTR  (i.e.,  with 
coolant  temperature  assumed  to  be  a constant)  following  Poore  (1973)  are 


dt 


C + Da(l-C)e^ 


(2.1  ) 


dT  T 

- T + BDa(1-C)e  - B (T-T  ) 
at  pc 


(2.2) 


where 


C = 


1-0,/C. 

A in 


T = 


„ (T-T.  ) 

E in 


B = 


rt: 


in 


(-AH)C.  E 
in 

RT^  pC 
in^  p 


Vk 


Da  = 


8 

p qpc. 


Here  the  dependent  variables  C,  T represent  dimensionless  conver- 
sion and  temperature.  The  parameters  (-AH),  k^,  E,  T^^,  Cj^,  R, 
Q»  Cp,  p,  V,  U,  A represent  heat  of  reaction,  rate  constant  of  reaction, 
activation  energy,  inlet  temperature  of  reactants,  reaction  concentra- 
tion, inlet  concentration  of  reactant,  universal  gas  constant,  reactant 
flow  rate,  reactant  specific  heat,  reactant  density,  volume  of  reactor, 
overall  heat  transfer  coefficient,  and  surface  area  for  heat  exchange 
respectively.  Here  we  have  made  the  exponential  approximation  (Aris, 
1975)  in  an  attempt  to  simplify  the  Arrhenius  kinetics.  Here  de- 
notes dimensioned  variables.  T^,  is  the  assumed  constant  dimensionless 
coolant  temperature  scaled  as  T.  If  we  neglect  this  assumption  and  take 


5 


6 


the  coolant  to  be  perfectly  mixed,  l.e.,  neglect  spatial  gradients,  the 
energy  balance  for  the  coolant  becomes 


dT 

Le  = a_(T  . -T  ) + (T-T  ) 
dt  2 cm  c c 


(2.3) 


where 


Le  = 


^c^c’^pc^  %^c^pc  „ 

’ “2  UA  ’ ^c 


E (T  -T.  ) 
c in 


UAV 


2 

RT7 

in 


T . 
cm 


E (T  . -T.  ) 
cm  m 


RT^ 

m 


The  derivation  of  these  equations  is  shown  in  Appendix  A.  A simi- 
lar energy  balance  for  the  cooling  jacket  was  done  by  Mukesh  and  Cooper 
(1986).  They  simulate  the  CSTR  using  the  outlet  coolant  temperature  as 
a variable  to  control  the  reactor.  Consequently  they  are  forced  to 
study  the  reactor  dynamics  by  including  the  energy  balance  for  the 
coolant.  They  study  the  effect  of  sampling  time  on  the  onset  of  limit 
cycles  in  their  reactor. 

In  (2.3)  Tq  is  the  outlet  coolant  temperature  which  is  the  same  as 
the  coolant  temperature  in  the  jacket.  However,  this  differs  from  the 
inlet  coolant  temperature  T^j^^  (see  Appendix  A).  We  have  introduced  two 
new  dimensionless  parameters  Le  and  a^.  The  parameter  Le  is  not  the 
same  as  the  traditional  Lewis  number.  It  is  similar  to  it  for  two 
reasons.  First,  it  multiplies  the  time  derivative  term  as  the  standard 
Lewis  number.  Secondly,  it  can  be  looked  upon  as  a ratio  of  two  time 
scales.  The  parameters  having  a subscript  'c'  represent  the  coolant 
variable  analogs  of  the  corresponding  reactor  variables.  The  equations 
modelling  the  CSTR  now  become  (2. 1-2.3). 
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The  choice  of  the  dimensionless  groups  is  such  that  we  have  a very 
close  parallel  to  Poore's  model.  The  parameters  B,  Da,  Sp  in  (2. 1,2. 2) 
are  similar  to  the  parameters  B,  Da,  B of  Poore's  model.  This  facili- 
tates comparison  of  our  results  with  Poore's  results  under  certain 
limiting  conditions  as  Le  ->■  0,  Le  ->■ 

The  advantages  in  choosing  the  dimensionless  groups  Le,  a in 
equation  (2.3)  are  as  follows: 

1.  The  parameter  Le  multiplies  the  time  derivative  term  in  (2.3). 
Hence,  the  steady-state  solution  obtained  by  setting  the  time 
derivative  to  zero  is  independent  of  Le.  This  facilitates 
computations  as  we  will  explain  later. 

2.  The  parameter  the  coolant  flow  rate  ocoure  in  only  one  dimen- 
sionless parameter  a^.  Hence,  can  be  used  as  a control 
parameter . 

3.  In  the  limits  of  Le  ^ 0,  Le  ^ the  results  offer  easy  compari- 
son with  Poore's  results.  Besides  in  these  limits  it  may  be 
possible  to  apply  singular  perturbation  methods  and  seek  the 
existence  of  relaxation  oscillations.  Our  system  may  also 
exhibit  quasi-per iodic  oscillations. 

The  steady-state  solutions  are  given  by 


+ ft  T 

0 = - C®®  + Da(1-C®®)e  ^ 


1 + 6 


M 


(2.4) 


T^®  = (BC^^ 


+ 6„T^.  )/(1 
M cm 


(2.5) 


nSS 


= (T 


ss 


" “2^ 


+ 


(2.6) 
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Here  superscript  'ss'  denotes  a steady  state  solution. 

" 1 ' + ct  is  ^ modified  or  effective  heat  transfer  coefficient.  The 
steady  state  solution,  thus,  depends  only  on  Da,  B,  It  is  indepen- 
dent of  Le  and  depends  on  6p,  through  6j^^.  Equation  (2.4)  is  identi- 
cal with  Poore's  equation  for  steady-state  conversion  with  his  6,  T 
playing  the  role  of  our  thus  use  Poore's  results  to 

map  his  steady  curves  to  our  curves  realizing  the  analogy  between  our 
3w  and  Poore's  6. 

The  nonlinear  equation  (2.4)  can  have  multiple  solutions  for  some 
combinations  of  B,  Da,  Following  Poore,  we  require  for  unique 

solutions  > 0.  Thus  a necessary  and  sufficient  condition  for  unique 

solutions  is  B ^ 4(1  + 6^^),  For  B >4(1  + B^^)  there  exists  a range  of 
Da  for  which  equation  (2.4)  has  three  solutions.  Figures  2.1  and  2.2 
are  plots  of  steady-state  conversion  versus  Da  for  the  two  cases  of 
unique  solutions  for  all  Da  and  multiple  solutions  for  some  values  of  Da 
respectively.  Typical  order  of  magnitude  estimates  for  the  various 
dimensionless  parameters  are  given  in  Table  1.  These  estimates  have 
been  obtained  from  the  values  of  parameters  quoted  in  experimental 
studies  in  the  literature  (see  Appendix  B).  All  computations  and 
diagrams  are  based  on  values  of  parameters  taken  from  these  ranges. 

Hence,  the  behavior  of  the  steady  states  is  very  similar  to  the 
standard  Poore  problem.  This  facilitates  the  study  of  the  role  of  the 
parameter  Le  on  the  onset  of  limit  cycles. 


(steady  conversion)  ^ C”  (steady  conversion) 
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2.1.  Unique  steady  states 


Figure  2.2.  Multiple  steady  states 
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TABLE  1 

Order  of  Magnitude 

Estimates  of  Dimensionless  Parameters 

Parameter 

Order  of  Magnitude 

B 

0 

1 

0 

0 

Da 

10'^  - 10“^ 

Le 

0 

1 

1 

0 

0 

“2 

0 

1 

1 

0 

0 

0 

0 

1 

1 

0 

CHAPTER  3 

LINEAR  STABILITY  ANALYSIS 


In  the  previous  chapter,  the  calculation  and  characterization  of 
the  steady-state  solutions  of  our  CSTR  were  described.  The  stability  of 
these  steady  states  to  infinitesimal  perturbations  is  now  discussed 
using  a linear  stability  analysis. 

Linearizing  the  equations  (2.1-3)  about  the  steady-state  solutions 
we  obtain 


X ' 


_d 

dt 


T' 

T' 

c 


-1 

cf 

1-C 

0 


- 6. 


-1  + BC®®  -6 


ss 


1-C 

Le 


ss 


(1  + OL^) 

Le 


x' 


T'  (3.1) 

T' 

c 


or 


_d 

dt 


with 


y'  = Cx^  T'.  T']^ 

where  x'  = - BC'  + T',  and  M is  the  Jacobian  matrix  occurring  in  (3.1). 

The  primed  variables  denote  perturbations  about  the  steady  state,  i.e., 

T =T  - T®®,  C'  = C - C®®,  T'  = Tq  - If  the  perturbations  grow 

c c 
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exponentially  with  time  having  a growth  constant  s,  i.e., 

^ ^ s t 

x'  = X e , we  obtain 

sy  = My  where  y = [x  , T , T^] 

The  stability  of  the  steady  state  is  now  determined  by  the  three 
eigenvalues  s of  the  matrix  M.  If  all  three  eigenvalues  have  negative 
real  parts,  i.e.,  lie  in  the  left  half  plane  of  the  complex  domain,  the 
steady  state  is  asymptotically  stable.  If  one  of  the  real  eigenvalues 
turns  from  negative  to  positive  as  we  vary  a parameter,  while  the  other 
two  eigenvalues  lie  in  the  left  half  plane  we  have  a steady-state  bifur- 
cation, i.e.,  the  solution  jumps  from  one  stable  steady  state  to  another 
steady  state.  If,  however,  the  real  part  of  a complex  conjugate  pair  of 
eigenvalues  changes  sign,  i.e.,  goes  from  negative  to  positive,  while 
the  other  eigenvalue  remains  negative,  the  steady  state  becomes  unstable 
and  we  see  the  onset  of  oscillations  or  limit  cycles.  The  dependent 
variables  now  vary  periodically  with  time.  Thus  a study  of  the  eigen- 
values 's'  of  M gives  us  information  about  the  stability  of  the  steady 
state  and  the  nature  of  the  onset  of  instability. 

The  eigenvalues  s of  M are  given  by  the  roots  of  the  cubic  charac- 
teristic equation 


3 2 

s - Trs  +Ps-D=0  (3.2) 

Here  Tr,  P,  D represent  the  trace,  sum  of  the  principal  (2x2) 
minors,  and  determinant  of  the  matrix  M.  The  expressions  for  Tr,  P,  D 


are 
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Tr 


(3.3) 


P 


B(1  + a^)  + 3 + (1  + a^) 


Le 


(3.4) 


i-Bz^  + Bz  - (1  + 6 )) 

M 


( 1 -z)Le 


(3.5) 


Here  the  variable  C®®  has  been  replaced  by  z for  convenience.  The 

expressions  for  Tr , P,  D are  all  quadratic  in  z.  The  expression  for  D 

is  exactly  like  Poore's  determinant  except  we  have  3 instead  of  3.  The 

M 

expression  for  Tr  is  identical  with  that  of  Poore's  trace  if  we  treat 


3 (with  3 playing  the  role  of  his  3).  Note  that  in  gen- 
eral  3.  Conditions  for  steady  state  bifurcation  and  Hopf 

bifurcation  will  now  be  given. 

The  necessary  and  sufficient  condition  for  asymptotic  stability  of 
the  steady  state  is  Re(s)  < 0 for  all  three  roots  of  equation  (3.2).  In 

terms  of  the  coefficients  of  the  cubic  this  translates  to  (Porter, 
1968), 


The  char ac ter  is t i c cubic  can  have  either  all  three  real  roots  or 
one  real  root  and  a pair  of  complex-conjugate  roots.  This  is  decided  by 
the  discriminant  of  the  cubic,  6,  where 


Tr  < 0,  P > 0,  D < 0,  TrP  - D < 0 
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6 = Up3  _ p^'Yr^  + 27D^  + 4DTr^  - iSTrPD 

For  6 < 0 we  have  all  three  real  roots  and  for  6 > 0 we  have  one 
real  root  and  a complex-conjugate  pair,  and  6=0  corresponds  to  three 
real  roots  of  which  at  least  two  are  equal. 

In  terms  of  the  cubic  coefficients,  the  conditions  for  steady-state 
bifurcation  is 


D = 0,  Tr  < 0,  P > 0.  (3.6) 
For  Hopf  bifurcation  the  corresponding  conditions  are 

TrP  - D = 0,  Tr  < 0,  P > 0,  and  6 > 0.  (3.7) 

Steady-state  bifurcation.  The  condition  for  steady-state  bifurca- 
tion D = 0 can  be  obtained  by  solving  the  quadratic  expression  in 
(3.5).  The  quadratic  is  similar  to  Poore's  expression  with  our  6 

M 

replacing  his  B-  Thus  for  a fixed  value  of  Bp,  ol^,  B,  we  calculate  the 
value  of  z and  hence  Da  at  which  this  bifurcation  occurs.  Note  that 
this  kind  of  instability  is  independent  of  Le.  In  order  to  observe  this 
instability  the  conditions  Tr  < 0,  P > 0 must  also  be  satisfied.  This 
ensures  we  have  two  eigenvalues  with  negative  real  parts  and  one  real 
eigenvalue  equal  to  zero. 

The  quadratic  expression  occurring  in  D is  identified  with  the 
quadratic  expression  arising  from  the  imposition  of  the  monotonic  varia- 
tion of  with  Da.  Thus  once  again  for  B < 4(1  + Bw)  we  have  no 
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steady-state  bifurcation.  For  B > 4(1  + the  quadratic  in  D has  two 
real  roots.  The  points  of  steady-state  bifurcation  hence  correspond  to 
the  turning  points  of  the  S-shaped  curve  of  C®®  versus  Da  provided 
conditions  (3.6)  hold.  The  occurrence  of  multiple  steady  states  gives 
rise  to  hysteresis  and  steady-state  bifurcation.  A sufficient  condition 
for  instability  arising  from  the  determinant  is  D > 0.  This  condition 
prevails  in  the  middle  branch  of  the  S-shaped  curve.  Note  that  this 
condition  is  the  opposite  of  the  condition  of  the  classical  CSTR  prob- 
lem, where  D < 0 is  sufficient  for  instability.  The  model  hence  has  a 
very  close  parallel  to  Poore's  model.  The  steady-state  behavior  does 
not  give  rise  to  any  new  phenomena. 

Hopf  bifurcation  and  the  onset  of  limit  cycles.  The  necessary  and 
sufficient  condition  for  Hopf  bifurcation  is  D < 0,  Tr  < 0,  P > 0, 

6 > 0 and  TrP-D  = 0.  The  condition  6 > 0 ensures  we  have  a complex 
conjugate  pair  of  eigenvalues,  and  TrP-D  = 0 ensures  their  real  part  is 
zero.  The  other  conditions  ensure  the  remaining  real  eigenvalue  is 
negative. 

The  condition  TrP-D  = 0 gives  rise  to  the  following  quartic  in  z 


-B(B  + BC1)z^  + ((B  + 1 + + C1)(B  + BC1 ) 

+ B(B  + BCl  + C2  + Cl)  - BCDz^  - (B(2C1  + 1+6  + C2) 

P 

+ (2  + Bp  + C1)(B  + BC1 ) + (B  + 1 ^ Bp  + C1)(B  + BC1  + Cl  + C2) 

(3.8) 

2 

- 2BC1)z  + ((B  + 1 + 6^  + C1)(2C1  +1+6  + C2)  + (2  + 6 + Cl ) 

P p p 
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(B  + BC1  + C2  + Cl)  - BC1  - (1  + B,JCl)z  - (2  + g + Cl  ) 

M p 


(2C1  + 1 + gp  + C2)  + (1  + gj^)C1 


where 


(1  + a ) g a 

Cl  = ^ ^ and  C2  = P 


Le 


Le 


This  can  be  rewritten  as  a quadratic  in  1/Le  as 


1 + a. 


(z-1)(Bz  - (B  + 1 + g^)z  + (2  + g^))( 


[(z-1)(Bz^  - Bz  + 1 + g ) + (Bz^  - (B  + 1 + g,Jz 

p M 


+ (2  + Bz  + (B  + 1 + gp)z  - (2  + gp)) 


(3.9) 


2 1 + a_ 

(-  Bz  + Bz  - (1  + g^))(1-z)]  {— -)  + 

M Le 


(-  Bz^  + (B  + 1 + g )z  - (2  + g ))(Bz^  - Bz  + 1 + g ) = 0 

P P P 


The  calculation  of  the  Hopf-bifurcation  points  using  the  quadratic 
(3.9)  is  mathematically  convenient.  It  is  made  feasible  because  the 
steady  state  is  independent  of  Le.  This  enables  us  to  fix  all  other 
parameters  in  the  problem  to  compute  the  critical  values  of  Le  at  which 
we  have  the  onset  of  instability.  The  volume  of  the  cooling  jacket 
occurs  only  in  Le  and  in  no  other  dimensionless  parameter.  Thus  Le  can 
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be  used  as  a control  parameter  to  study  the  dynamics  of  the  CSTR. 
Varying  Le  would  be  equivalent  to  changing  the  coolant  residence  time. 

Ray  and  Hastings  (1980)  modelled  the  CSTR  by  considering  the  ther- 
mal capacitance  and  material  capacitance  of  the  reactor  to  be  different 
from  those  of  the  reacting  fluid.  This  more  general  model  contains  a 
similar  parameter  'Le'  in  their  energy  balance.  They  find  that  their 
model  has  periodic  solutions  for  sufficiently  low  values  of  Le  their 
Lewis  number.  Their  parameter  Le  is  different  from  our  definition  of 
Le.  It  however  arises  because  of  considering  an  additional  thermal 
capacitance.  They  exploit  successfully  the  independence  of  the  steady 
state  of  the  Le  number  to  1 ) calculate  critical  values  of  Le  for  Hopf 
bifurcations  and  2)  study  relaxation  oscillations  under  the  limit  of 
Le  -*■  0. 

The  two  limiting  cases  of  interest  are  Le  -»•  “ and  Le  ->•  0.  These 
allow  an  easy  comparison  of  our  results  with  Poore's  results.  Under 
these  asymptotic  conditions  we  derive  conditions  for  the  Hopf-bifurca- 
tion  points. 

Case  of  Le  ^ °°.  Under  these  conditions  we  divide  equation  (2.3)  by 
dT 

Le,  to  get  = 0 or  T^  = T^(t=0),  a constant.  Here  the  thermal  iner- 
tia of  the  cooling  jacket  is  sufficiently  high  and  prevents  the  coolant 
temperature  from  rising  significantly  past  its  initial  value.  Thus 
equation  (2.3)  is  decoupled  from  (2.2)  and  T^  can  be  regarded  a constant 
TQ(t=0)  in  equation  (2.2). 

The  Hopf  bifurcation  points  can  then  be  obtained  for  this  case  from 
the  standard  Poore  model  if  we  use  our  6p  for  Poore's  6.  We  obtain  the 
same  result  by  taking  the  limit  Le  ->•  ® in  the  quartic  to  obtain  the 


quadratic. 
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2 

Bz  -(B+l+B)z+(2+B)=0 
P p 

Case  Le  -»•  0.  Under  this  limit  we  make  the  pseudo-steady  state 
approximation,  i.e.,  we  take  equation  (2.3)  to  be  always  at  steady 
state.  Under  these  conditions  becomes  a fast  variable  (Mishchenko 

and  Rozov,  1980).  We  eliminate  between  (2.2)  and  (2.3).  The  dynam- 

ics of  the  reactor  are  then  identical  with  those  of  the  standard  Poore 
reactor  with  his  6 now  being  replaced  by  B^.  Now  the  Hopf  bifurcation 
points  occur  at  the  roots  of  the  quadratic, 

Bz^  - (B  + 1 + R)z  + (2  + B„)  = 0 
M M 

This  is  identical  with  Poore's  Hopf  bifurcation  condition  with  B 
replaced  by  Bj^.  We  thus  see  that  under  the  two  limits  of  Le  ^ 0,  Le  ->  «> 
we  recover  the  standard  Poore  problem  for  different  values  of  B. 

The  effect  of  the  finite  thermal  capacity  (Le)  of  the  cooling 

jacket  on  the  asymptotic  stability  of  the  steady  states  can  be  described 

using  plots  of  LOq  versus  Da.  Here  Loq  represents  critical  values  of 

Le,  i.e.,  Le^  is  the  root  of  (3.9),  when  we  fix  z,  a„,  B,  B • 

2 P 

The  quadratic  (3.9)  in  Le  can  have  a maximum  of  two  positive  real 
roots  LOq^  , LOq2  for  some  Da.  Let  LeQ-]  < LeQ2.  For  values  of  Le,  Le  < 
LSqp  Le  > Le^2  we  have  unstable  steady  states  for  those  Da.  For  Le^^  < 
Le  < LeQ2>  we  have  asymptotically  stable  steady  states  for  these  Da. 
Instability  under  the  conditions  of  Le  0,  Le  ® is  to  be  expected 
since  these  asymptotic  limits  correspond  to  Poore's  problem  under  the 
limit  of  B = B = Bp  respectively.  We  can  verify  these  results  from 


19 


Figure  3 


1.  Critical  Lewis  numbers  for  stable  and  unstable  states 
of  conversion  shown  in  Fig.  3.2 


Figure  3.2.  Steady  state  conversion 
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Figure  3.3.  Critical  Lewis  numbers  for  stable  and  unstable 
states  of  conversion  shown _in  Fig.  3. A 


Figure  3.4.  Steady  state  conversion 
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the  stability  of  Poore’s  model  for  these  values  of  3.  Figure  3.2  depicts 
the  steady-state  conversion  as  a function  of  Da  for  a fixed  B, 

^p’  ^^Sure  3-1  is  the  plot  of  the  critical  Le  versus  Da  for  the 
steady  states  shown  in  Fig.  3*2.  The  lower  curve  corresponds  to  Lep. , 
the  upper  curve  to  LeQ2»  values  of  Da,  where  teQ^  exists  the  steady 

states  are  unstable  for  Le  < LeQi  . If  we  design  our  reactor  with  Le  > 
Le^2  the  steady  state  is  again  unstable.  Choosing  Le  such  that  LC(^.]  < 
Le  < LeQ2  for  the  Da,  B,  6^,  shown  assures  us  of  an  asymptotically 
stable  steady  state.  Thus  by  appropriately  choosing  Le  we  can  operate 
the  reactor  at  an  asymptotically  stable  steady  state. 

In  Fig.  3.3  is  a similar  plot  of  Le^  versus  Da  for  a fixed 
B,  3p,  a^.  The  corresponding  steady-state  conversion  is  shown  in  Fig. 
3.^.  Here  the  combination  of  parameters  is  such  that  the  equation  (3.9) 
has  only  one  positive  real  root  LOq.  For  Le  > LeQ  steady  state  for  the 
corresponding  Da  shown  in  Fig.  3-^  is  asymptotically  stable.  At  Le  < 
Le^  we  have  unstable  steady  states.  Under  such  cases,  we  can  operate  at 
stable  steady  states  by  having  sufficiently  high  values  of  Le.  Ray  and 
Hastings  (1980)  have  a similar  result.  They  find  onset  of  limit  cycles 
for  sufficiently  low  values  of  their  Le.  Their  definition  of  Le  is 
quantitatively  different  from  ours.  It  represents  a measure  of  thermal 
inertia  of  the  reactor  and  is  hence  similar  physically. 


CHAPTER  4 

THE  D-PARTITION  METHOD  AND  APPLICATIONS  TO  THE  CSTR  MODEL 

The  first  order  exothermic  reaction  has  been  investigated  by  many 
researchers  in  great  detail.  Most  publications  (Poore  (1973),  Uppal  et 
al . (1974,  1976))  to  date  have  been  concerned  with  the  evaluation  of  the 
stability  of  a steady  state  for  a fixed  combination  of  parameter 
values.  This  approach  to  the  problem,  however,  has  some  limitations. 
This  problem  is  highly  nonlinear  and  requires  detailed  calculations. 
Moreover,  the  resultant  steady-state  conversion  of  the  reactor  arises  as 
a result  of  the  problem  and  is  not  predetermined.  In  this  section  we 
wish  to  address  the  inverse  problem,  i.e.;  what  are  the  total  possible 
combinations  of  parameter  values  which  give  us  a desired  level  of  con- 
version? Which  of  these  parameter  combinations  yield  an  asymptotically 
stable  fixed  level  of  conversion? 

In  an  industry  it  is  usually  of  interest  to  run  a reactor  at  a 
fixed  level  of  conversion.  It  is  hence  necessary  to  see  what  combina- 
tion of  parameter  values  yield  this  particular  conversion.  We  can  then 
design  our  reactor  appropriately.  An  elegant  way  of  doing  this  is  by 
using  the  D-partition  method  (Porter,  1968).  We  first  describe  the 
method,  its  salient  features  and  its  limitations.  We  conclude  the 
chapter  by  applying  it  to  a couple  of  situations. 

The  D-partition  method  can  be  used  effectively  to  deal  with  charac- 
teristic equations  whose  coefficients  depend  on  certain  parameters  p 
(P-|.  •••»  Pn^  linearly.  Here  p^  , P2»  etc.  are  scalars.  The  method 
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is  used  to  divide  the  p parameter  space  into  different  regions.  It 
gives  us  information  on  how  the  roots  of  the  characteristic  equation 
change  as  we  go  from  one  region  to  another.  In  particular,  we  can 
identify  the  number  of  roots  with  negative  real  parts  in  each  region. 
We  then  obtain  the  combination  of  parameter  values  where  all  roots  of 
the  characteristic  equation  have  negative  real  parts  Indicating  asympto- 
tic stability.  The  application  is  particularly  elegant  to  characteris- 
tic equations  dependent  on  two  parameters,  as  we  can  now  depict  graphi- 
cally the  different  regions  in  a two  dimensional  diagram. 

Let  the  characteristic  equation  be 

F(s)  = 0 (4.1) 

where  F(s)  is  usually  an  nth  order  polynomial  in  s.  Here  s is  a complex 
variable  (o  + iw).  Over  the  complex  domain  the  region  of  asymptotic 
stability  is  the  left-half  plane,  the  region  where  all  roots  have  nega- 
tive real  parts.  The  D-partition  method  maps  the  boundary  between  the 
regions  of  stability  and  instability  in  the  s-plane,  i.e.,  the  entire 
imaginary  axis  (s  = iw)  into  the  parameter  space  of  interest.  This 
divides  the  entire  parameter  space  into  different  regions.  Each  region 
can  be  denoted  by  D(k,n-k),  0 ^ k $ n,  where  k (n-k)  represents  the 
number  of  roots  of  equation  (4.1)  with  negative  (positive)  real  parts, 
and  n the  order  of  polynomial.  The  region  of  asymptotic  stability  hence 
corresponds  to  the  region  D(n,0),  where  we  have  no  roots  with  positive 
real  parts.  If  we  neglect  the  possibility  of  a characteristic  root 
Jumping  from  to  +“  an  increase  in  number  of  roots  with  positive  real 
parts  can  occur  only  as  a result  of  a root  crossing  the  imaginary  axis 
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from  the  left  to  the  right  half  plane.  A root  crossing  the  imaginary 
axis  from  the  left  to  the  right  half  plane  as  we  vary  parameters  corre- 
sponds to  a representative  point  in  the  parameter  space  moving  from 
region  D(k,n-k)  to  region  D(k-1 , n-k+1).  Thus  for  finite  roots  the  D 
partition  boundaries  are  the  maps  of  the  imaginary  axis  (s  = iw), 

-00  < w < +00  in  parameter  space. 

This  method  is  efficient  and  informative  in  constructing  domains  of 
stability  in  terms  of  parameters  occurring  in  the  problem.  Although  it 
is  possible  to  construct  these  domains  in  more  than  two  dimensional 
spaces  it  can  be  applied  elegantly  to  problems  containing  one  or  two 
parameters.  The  method  of  construction  of  stability  domains  for  the  one 
and  two  parameter  cases  will  now  be  demonstrated. 

— £^-^meter case . If  the  coefficient  of  the  characteristic 

equation  were  to  depend  on  one  parameter,  say  5,  linearly,  we  can 
rewrite  equation  4.1  as 


P(s)  + 5Q(s)  = 0 


or 


The  values  of  C corresponding  to  the  imaginary  axis  in  the  s-plane  is 
obtained  by  substituting  s=iw  in  (4.2) 


5(w) 


P(iw) 

Q(iw) 


-“  < w < ® 


(4.3) 
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This  equation  now  defines  the  D-partition  boundary  corresponding  to 
finite  roots.  We  let  w vary  from  to  and  follow  the  locus  traced 
by  equation  (4.3),  in  the  Re(5)  versus  Im(C)  plot.  If  we  shade  the  left 
hand  side  of  this  locus  as  w increases  from  -=»,  then  the  shaded  side 
corresponds  to  the  more  stable  side  of  the  boundary.  Thus  crossing  a 
boundary  from  the  shaded  to  the  unshaded  side  corresponds  to  a root 
passing  from  the  left  to  the  right  half  of  s-plane. 

If  all  the  boundaries  are  shaded  according  to  the  above  convention, 
we  can  determine  the  value  of  k for  each  domain  D(k,n-k)  once  we  know 
the  value  of  k for  a domain.  This  is  done  by  observing  the  number  of 
times  a boundary  is  crossed  and  decreasing  or  increasing  the  value  of  k 
according  as  the  crossing  is  from  the  shaded  to  the  unshaded  side  or 
vice  versa.  We  thus  obtain  the  region  of  asymptotic  stability  D(n,0). 
For  a real  parameter  C,  this  asymptotically  stable  region  reduces  to  the 
segments  of  the  real  axis  in  D(n,0). 

The  two  parameter  case^  We  now  consider  the  case  where  the  coeffi- 
cients of  the  characteristic  equation  depend  linearly  on  two  parameters 
n.  The  equation  (4.1)  can  now  be  rewritten  as 


F(s)  = CP(s)  + nQ(s)  + R(s)  = 0 


(4.4) 


Taking 


P(iw)  = P^(w)  + iP2(w) 
Q(iw)  = (w)  + 


26 


R(iw)  = (w)  + iR^Cw) 

Along  the  imaginary  axis  (s=iw),  we  obtain 

F(iw)  = 5P^(w)  + nQ^(w)  + R^(w)  + i(cP2(w)  + nQ^Cw)  + R2(w))  = 0 
This  implies 


(w)  + nQ^ (w)  + R^ (w)  = 0 


^P^Cw)  + nQ2(w)  + R2(w)  = 0 


(4.5) 


The  D-partition  boundaries  corresponding  to  finite  roots  for 


A(w)  = P^Q^  - P^Q^  ^ 0 


(4.6) 


is 


5 = (-  R^Q2  + R2Q^)/A 

n = (-  R2P^  + R^P2)/a. 


(4.7) 


For  A(w)  * 0 As  w increases  from  -»  to  +®  these  equations  trace  out 
a locus  in  the  plane  which  corresponds  to  a map  of  the  imaginary 
axis  in  the  s-plane.  This  is  called  the  non-singular  locus  and  is  a 
part  of  the  system  of  D-partition  boundaries. 
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(a) 
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Figure  4:1.  Shading  rules  for  D-partition  boundaries.  Locii 
intersecting  at  (a)  v ^ 0;  (b)  w = 0 
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If  A(w)  = 0 for  some  w , we  obtain  the  corresponding  singular  locus 
by  substituting  the  value  of  w in  (^.5). 

The  various  domains  D(k,n-k)  0 $ k ^ n must  now  be  identified.  To 

do  this  we  need  a convention  to  shade  the  various  boundaries.  If  the 

axis  in  the  plane  constitutes  a right-handed  orthogonal  system,  the 
non-singular  locus  must  be  shaded  on  the  left  for  A > 0 and  on  the  right 
for  A < 0 if  the  locus  is  traversed  in  the  direction  of  increasing  w. 
It  is  usually  found  that  the  same  curve  is  traversed  twice,  when  w 
increases  from  -“  to  0 and  when  w increases  from  0 to  ” and  A reverses 
sign  at  w = 0,  ± ®.  Here  the  non-singular  locus  will  be  shaded  twice  on 
the  same  side.  Crossing  such  a double  shaded  boundary  corresponds  to 

two  characteristic  roots  crossing  the  imaginary  axis  simultaneously  as 
we  vary  the  parameters. 

The  rules  for  shading  a singular  locus  depend  on  the  value  of  w at 
which  it  intersects  the  non-singular  locus  (Porter,  1968).  If  the  two 
locii  intersect  at  a finite  non-zero  value  of  w,  at  which  A changes 
sign,  the  shading  is  done  as  in  Fig.  4.1a.  If  the  intersection  occurs 
at  a point  on  the  non-singular  locus  at  which  A changes  sign,  at  w = 
0,  ± “,  the  shading  is  done  as  shown  in  Fig.  4.1b.  Here  the  shading  is 

on  the  same  side  for  both  locii  hear  the  point  of  Intersection,  but 

beyond  this  point  the  shading  passes  to  the  other  side  of  the  singular 
line.  If  A vanishes  identically,  the  singular  lines  are  the  only  bound- 
aries. 

Applications . In  this  section  we  discuss  the  application  of  the  D- 
partition  method  to  two  chemical  engineering  problems. 

(1)  Classical  CSTR.  The  classical  CSTR  system  as  modelled  by  Poore 
(1973)  with  a constant  coolant  temperature  has  a characteristic  equation 


given  by 


29 


2 2 

F(s)  = s + (Bz  - (B+1+6)z  + (2+B)) 


(Bz^  - Bz  ( 1 +3) ) 
( 1 -z) 


0 


(^.8) 


Here  z represents  conversion.  B and  3 are  the  dimensionless  param- 
eters representing  heat  of  reaction  and  heat  transfer  coefficient  as 
defined  by  Poore.  For  a fixed  z the  steady  state  conversion  level,  F is 
now  linear  in  B and  3.  Hence  the  method  of  D-partition  can  be  applied 
to  calculate  stability  domains  elegantly  in  (B— 3)  space. 

Using  the  method  explained  in  the  previous  section,  we  find  the 
singular  locus  corresponds  to  w=0  and  is 


The  non-singular  locus  corresponds  to 


z 


2 

Eliminating  w , this  reduces  to 


3 > 3 


c 


z 


(4.10) 


B > B 
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The  two  locii  intersect  at  6^,  and  have  been  drawn  in  Fig.  4.2  for  z 
= .5  and  shaded  accordingly.  The  two  locii  (4.9,  4.10)  are  straight 
lines.  The  equation  for  the  singular  and  non-singular  locus  is  the  same 
as  the  equation  obtained  by  setting  the  determinant  and  trace  of  the 
Jacobian  matrix  to  zero  for  a fixed  z.  They  divide  the  B-S  parameter 
space  into  three  regions  A,B,C.  Region  A corresponds  to  B-B  combina- 
tions where  the  characteristic  equation  has  two  eigenvalues  with  nega- 
tive real  parts.  In  region  B,  one  eigenvalue  has  a positive  real  part 
and  the  other  a negative  real  part.  In  region  C both  characteristic 
roots  have  positive  real  parts.  Thus  B-6  values  in  region  A give  rise 
to  a reactor  which  can  operate  at  a fixed  stable  steady  state  conversion 
level  of  50%  if  Da  is  chosen  according  to 

Bz 

Da  = ^ ^ K (4.11) 

We  have  used  the  D-partition  method  to  divide  our  parameter  space 
into  different  regions.  The  method  is  elegant  since  the  three  dimen- 
sionless parameters  Da,  B,  6 in  the  problem  can  be  determined  explicitly 
for  a fixed  stable  conversion  level.  The  stability  domain  has  been 
drawn  in  the  B-g  parameter  space.  No  nonlinear  equations  need  to  be 
solved  in  this  method  since  we  can  use  (4.11)  to  obtain  Da  explicitly 
for  a fixed  z,  B,  6.  The  design  of  a CSTR  for  a fixed  conversion  level 
can  hence  be  done  elegantly  using  this  method. 

Our  ability  to  eliminate  the  steady-state  temperature  in  favor  of 
conversion  using  a linear  transformation  makes  the  method  powerful. 
This  enables  us  to  fix  the  conversion  of  our  reactor  at  the  desired 
level.  The  temperature  is  now  fixed  automatically  at  a given  level  and 
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13 


Figure  4.2.  D-partition  diagram  for  the  classical  CSTR, 
Poore  model  (1973) 
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we  do  not  have  to  concern  ourselves  with  it  as  long  as  it  is  not  unreal- 
istic. If  we  were  unable  to  eliminate  the  steady-state  temperature 
using  our  linear  relationship,  we  would  still  be  able  to  use  the  above 
method.  The  method,  however,  would  only  be  locally  valid  for  a certain 
fixed  temperature  and  conversion.  Our  ability  to  eliminate  temperature 
in  favor  of  conversion  makes  our  analysis  globally  valid. 

(2)  CSTR  with  coolant  temperature  non-constant.  Here  for  our 
model  equation  (2.1-3)  the  character istic  equation  (3.1)  is  a cubic. 


3 2 

s - Trs  + Ps  - D = 0 


(4.12) 


The  coefficients  Tr , P,  D contain  the  parameters  z,  a , B,  6 , Le  (3.3- 

P 

3.5).  These  parameters  occur  in  nonlinear  combinations  with  each 

other.  It  appears  hence  that  the  D-partition  method  is  inapplicable. 

However,  by  fixing  some  of  the  parameters  at  fixed  values,  we  can  treat 
4.12  to  be  linear  in  the  remaining  one  or  two  parameters.  This  enables 
us  to  determine  stability  domains  in  terms  of  these  parameters  of  inter- 
est . 

In  (4.12),  by  assigning  values  to  z,  a^,  Le,  we  can  treat  the  equa- 
tion to  be  linear  in  B,  6p.  We  can  then  evaluate  stability  domain 
in  B-g  space. 

Fixing  z,  a^,  B,  0^  enables  us  to  treat  the  cubic  equation  as 
linear  in  1/Le.  Stability  domains  in  the  Re  (7“]  versus  Im( — -)  plots 

Ljc  L0 

can  then  be  obtained.  With  a view  to  verifying  the  Figs.  3.1,  3.3  the 
D-partition  diagrams  for  (4.12)  were  drawn  in  terms  of  1/Le.  Figure  4.3 
is  the  D-partltion  diagram  which  corresponds  to  Fig.  3.1.  Here  B, 
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Figure  4.3.  D-partition  diagram  in  1/Le  plane  of 
equation  (4.12) 
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-2.0 


Figure  4.4.  D-partition  diagram  in  (1/Le)  plane  of 
equation  (4.12) 
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6p»  0^2  been  fixed  at  the  same  values  as  in  Fig.  3*1  « We  choose  a z 

and  hence  a Da  from  Fig.  3.2,  at  which  we  have  a Le^^  , LeQ2  in  Fig.  3.1 

and  obtain  the  corresponding  region  of  stability  A,  D(3,0).  Since  the 
parameter  Le  is  real,  the  values  of  Le  ensuring  stability  for  the  chosen 

steady  state  is  the  section  of  the  real  axis  in  D(3,0).  This  is  bounded 

on  the  real  axis  between  Le^^  , LeQ2*  Thus  the  region  of  stability  is 
Lsq^  < Le  < LeQ2  as  found  in  Fig.  3.1. 

For  the  D-partition  diagram  in  Fig.  the  parameters  B,  $ , 

were  chosen  to  be  the  same  as  those  in  Fig.  3.3.  Choosing  a z from  Fig. 
3.4,  we  find  the  corresponding  region  of  stability  as  the  intercept  of 
D(3i0)  on  the  real  axis,  i.e.,  Le  > LSq.  The  other  intercept  of  the 
boundary  passes  through  the  origin  or  Le  = ®>.  This  confirms  the  result 

of  Fig.  3.3. 

The  D-partition  method  can  be  used  to  design  a cooling  Jacket  for 
the  stable  operation  of  a reactor  at  a given  steady  state.  This  is 
possible  since  the  parameter  occurs  only  in  Le  and  in  no  other  para- 


meter . 


CHAPTER  5 

CONCLUSIONS  AND  SCOPE 


The  use  of  different  nonlinear  models  in  chemical  engineering  is  an 
old  concept.  In  this  dissertation  the  standard  model  for  the  CSTR  has 
been  modified  and  extended  to  a more  realistic  model  by  considering  the 
finite  thermal  capacity  of  the  cooling  jacket.  The  D-partition  method 
and  its  applications  to  the  standard  CSTR  model  and  the  model  developed 
here  have  been  studied.  The  different  phenomena,  like  multiple  steady 
states,  limit  cycles  observed  in  Poore's  model  of  the  CSTR  occur  over  a 
very  narrow  range  of  parameter  values.  It,  hence,  appears  that  these 
phenomena  are  more  an  exception  than  the  rule,  and  the  natural  charac- 
teristic of  the  reactor  is  to  have  stable  steady  states.  The  finite 
thermal  capacity  of  the  cooling  jacket  was  included  in  the  model  to  see 
if  the  region  in  parameter  space  exhibiting  limit  cycles  became  wide  or 
not.  It  appears  that  the  region  over  which  we  see  Hopf  bifurcation  does 
decrease  and  hence  the  finite  thermal  capacity  has  a stabilizing  influ- 
ence on  the  steady  states. 

The  disappearance  of  the  Hopf-bifurcation  points  for  some  values  of 
Le  raises  the  question — do  we  have  an  isola  of  periodic  solutions?  And 
if  we  do,  are  there  any  methods  of  predicting  such  isolas?  Preliminary 
simulations  of  our  equations  show  the  presence  of  quasi-per iodic  solu- 
tions to  our  equations  (2. 1-2.3). 

The  parameter  Le  multiplies  a time  derivative  term.  Under  the 
limits  of  Le  > 0,  Le  > ®,  we  arrive  at  singular  perturbation  limits. 
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Under  these  limiting  conditions,  we  have  two  vastly  differing  time 
scales.  Some  of  the  equations  are  always  at  steady  state.  Here  some  of 
the  variables  experience  rapid,  almost  instantaneous,  variations.  The 
periodic  solutions  are  discontinuous  and  now  exhibit  successive  alterna- 
tion between  fast  and  slow  motion.  This  kind  of  motion  is  called  a 
relaxation  oscillation.  It  would  be  important  to  study  the  dynamics  of 
our  system  under  different  limits  of  parameter  values  to  check  for  such 
oscillations.  Such  oscillations  have  been  observed  in  many  systems 
(Aluko  and  Chang,  1986). 

In  the  second  part  of  the  dissertation,  the  D-partition  method  has 
been  described.  This  method  can  be  used  to  calculate  stability  domains 
in  terms  of  the  different  parameters  occurring  in  the  model.  Most 
models  contain  a large  number  of  parameters.  These  usually  occur  in 
nonlinear  combinations  with  each  other.  In  the  mathematical  analysis  of 
these  models  the  effect  of  a parameter  on  the  system  dynamics  is 
studied.  The  parameter  is  chosen  usually  from  the  point  of  view  of 
mathematical  convenience  and  is  usually  not  a practical  parameter  which 
can  be  varied.  In  the  D-partition  method,  we  can  study  and  choose  a 
practically  convenient  parameter  for  the  design  of  any  system.  This  can 
be  achieved  by  fixing  the  different  parameters  in  the  model  at  reason- 
able assigned  values  and  constructing  stability  domains  about  the  one  or 
two  parameters  of  interest  which  are  chosen  as  control  parameters. 

Extension  to  flickering.  There  have  been  numerous  studies  on  the 
onset  of  turbulence  in  natural  convection  (Howard  and  Kr ishnamurthy , 
1986).  These  studies  have  focused  on  the  different  routes  in  the  tran- 
sition to  turbulence.  In  natural  convection,  close  to  initiation  the 
induced  fluid  motion  is  very  regular.  Under  very  high  temperature 
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gradients  the  ordered  fluid  motion  breaks  down  and  gives  rise  to  chaotic 
motion  in  the  fluid  phase.  This  chaotic  motion  induces  a chaotic  varia- 
tion of  the  surface  temperature.  The  phenomenon  of  "flickering"  or 
random  temperature  fluctuations  has  been  observed  on  catalytic  wires  or 
gauzes.  It  is  possible  that  these  are  induced  by  natural  convection. 
The  initiation  of  natural  convection  induced  by  the  exothermicity  of  a 
surface  reaction  has  been  studied  by  Viljoen  and  Hlavacek,  1987.  We 
have  to  resort  to  complex  numerical  schemes  to  study  the  dynamics  of 
these  partial  differential  equations.  The  next  step  would  be  to  retain 
the  essential  features  of  the  model  and  simplify  the  mathematical  equa- 
tions to  get  a tractable  system  of  equations. 

At  the  initiation  of  natural  convection  cells  are  induced  on  the 
surface  (Chandrasekhar,  1961).  This  induces  a regular  pattern  of  alter- 
nating hot  spots  and  cold  spots  on  the  catalytic  surfaces.  These  are  in 
thermal  communication  with  each  other  through  the  fluid  and  thermal 
conductivity  of  the  catalyst.  A simplified  representation  of  this  model 
would  be  a network  of  CSTRs  each  communicating  with  each  other.  This  of 
course  has  also  the  advantage  of  being  a lumped  model.  The  first  step 
towards  studying  such  a network  would  now  be  to  examine  two  CSTRs  where 
the  effluents  from  one  is  used  as  a coolant  in  the  other  and  vice 
versa.  This  system  has  four  dependent  variables. 

It  has  been  known  for  a long  time  (Lorenz,  1963)  that  a system  of 
three  nonlinear  ordinary  differential  equations  although  completely 
deterministic  gives  rise  to  chaotic  motions.  Since  then  there  have  been 
numerous  models  in  a variety  of  disciplines  exhibiting  chaos.  This  has 
led  people  to  believe  that  chaotic  dynamics  for  a nonlinear  system  is 
more  the  rule  than  the  exception.  With  a view  to  studying  a system  with 
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the  minimum  three  dependent  variables,  the  two  CSTR  networks  described 
above  can  be  further  simplified.  The  coolant  temperature  can  be  used  as 
a third  state  variable  in  the  modelling  of  a single  CSTR.  The  similar- 
ity between  equations  governing  the  CSTR  and  those  modelling  a surface 
reaction  has  been  known  for  a long  time  (Aris,  1975).  Under  the  limit 
of  high  thermal  conductivity  and  low  heat  transfer  coefficient,  the 
equations  become  identical,  as  we  can  now  neglect  spatial  variations  and 
the  two  systems  are  described  as  lumped  systems  (see  Appendix  C).  If  we 
take  into  account  the  heat  gained  by  the  reactant  fluid  due  to  the 
exothermici ty  of  the  surface  reaction,  the  modelling  equations  of  the 
surface  reaction  become  identical  to  equations  (2.1-3).  Thus  any  cha- 
otic dynamics  exhibited  by  our  system  can  be  possibly  used  to  interpret 
the  occurrence  of  flickering.  It  would  thus  be  of  interest  to  study  the 
possibility  of  equations  (2.1-3)  exhibiting  chaotic  motions. 


APPENDIX  A 

DERIVATION  OF  THE  GOVERNING  LUMPED  EQUATIONS  OF  A CSTR 


In  this  appendix  the  governing  dimensionless  equations  for  the  CSTR 
are  derived.  In  general  the  coolant  temperature  varies  along  the  cool- 
ing coil.  The  lumped  energy  balance  for  the  coolant  is  derived  rigor- 
ously using  a formal  perturbation  argument.  The  reactor  is  assumed  to 
be  well  mixed.  We  neglect  spatial  variations  in  the  reactor. 

The  mass  balance  equation  written  for  the  reactor  is 


V 


dt 


-E/RT 


(A.1  ) 


The  energy  balance  equation  written  for  the  reactor  is 


^ " (-AH)VkC^e'^^^^  - UA(T-T^) 


(A. 2) 


We  represent  the  cooling  coil  as  a tubular  reactor  with  no  reaction 
occurring  in  it.  The  energy  balance  for  the  coil  is  then 


3T  3T 

-sf  ■ ><1  ^ ^ 


subject  to 


3T^(0) 

3x 


e(T  (0)  - T 
c 


cm 


9T^(1  ) 
9x 


0 


(A. 3) 
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where 


e 


(v  p C 
c c p 


T 


All  the  variables  are  explained  in  the  Notation.  Here  T varies  with 

c 

space  and  time.  The  space  coordinate  'x'  along  which  T varies  is 

c 

measured  along  the  cooling  coil.  We  seek  the  form  of  these  equations 
under  the  limit  e ^ 0 i.e.  infinite  thermal  conductivity  of  the  fluid. 
We  expand  the  in  a formal  perturbation  series  as 


T = Z T .e' 
" i=0 


(A. 4) 


Substituting  (A. 4)  into  (A. 3)  and  equating  powers  of  t,  we  obtain 


0(e°)  ^ = 0 with 

9x 


9T^^(0)  9T  (0) 

_ CO  = 0 = 


3x 


ax 


(A. 5) 


This  yields  T^^  is  independent  of  x.  Its  variation  with  time  is 


determined  by 


2~ 

9T  9 T 

‘■oS  -5T  ■ — r - 

c ax 


(A. 6) 


ax 


= T (0)  - T . 
CO  cm 


aT^. (1 ) 

Cl 

ax 


0 
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Integrating  (A. 6)  with  respect  to  x and  using  the  boundary  condi- 
tions and  remembering  is  independent  of  x,  we  arrive  at 


dT 

CO  — ^ - 

p C V — : — = (T  -T  ) - UA(T-T  ) 
^c  p^  c dt  ^ cin  co'^  co"^ 


(A. 7) 


To  a first  order  approximation,  i.e.,  in  the  limit  e ^ 0,  we  can 


approximate  the  coolant  jacket  temperature  as  T (x,t)  = T (t) 

C CO 


Defining  a dimensionless  temperature  T as  T = — (T-T  ) we  can 

,2  in 


rewrite 


RT. 

in 


RT. 

-E/RT.  1 + T 
^-E/RT  in  ^ E 

e = e .e 


-E/RT. 

in  T 
= e .e  . 


This  is  the  exponential  approximation.  It  is  accurate  in  the 

limit  E ^ Under  this  limit  1 « T(RT.  )/E.  If  we  now  define  the 

in 

dimensionless  parameters  B,  a^.  Le,  Da  as  in  the  text,  equations 
(A.1,  A. 2,  A. 7)  reduce  to  equations  (2. 1-2. 3). 


appendix  b 

VERIFICATION  OF  THE  MODEL  WITH  EXPERIMENTS  IN  THE  LITERATURE 


There  have  been  a few  experimental  studies  on  the  onset  of  oscilla- 
tions. Baccaro  et  al.  (1970)  studied  the  hydrolysis  of  acetyl  chloride, 
Bush  (1969)  studied  the  vapor  phase  chlorination  of  methyl  chloride  and 
Vermeulen  et  al.  (1986)  studied  the  hydration  of  2,3  epoxy-1  propanol 
catalyzed  by  sulfuric  acid.  In  this  appendix  the  experimental  data 
presented  in  the  above  mentioned  papers  have  been  used  to  compute  the 
dimensionless  parameters  occurring  in  our  -model.  These  values  are 
presented  in  Table  B1  . The  order  of  magnitude  estimates  for  the  various 
dimensionless  parameters  in  Table  1 of  the  text  are  based  on  these 
values . 


Table  B1  . Dimensionless  parameters  calculated  from  experimental 


data  in  literature 

Baccaro  et  al. 
(1970) 

Bush 

(1969) 

Vermeulen  et 
(1986) 

al. 

B 

14.8 

40.0 

31  .7 

6 

P 

5.3 

2.0 

6.3 

Da 

.14 

— 

.02 

Le 

.06 

— 

— 

“2 

.77 

— 

— 

data  of  Baccaro 

et  al 

. is  complete 

as  we  know 

the  values  for  all 

parameters  from 

their 

experiments. 

For  the  set 

of  parameters 

of 
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Baccaro  et  al.  our  model  predicts  oscillations  around  a mean  value  of 
17°C.  They  observe  oscillations  around  22°C. 


APPENDIX  C 

EQUIVALENCE  OF  THE  EQUATIONS  GOVERNING  A CSTR  AND 
THOSE  OF  A SURFACE  REACTION  ON  A WIRE 


In  this  appendix  we  establish  the  similarity  between  the  equations 
modelling  our  CSTR  and  those  modelling  a surface  reaction  on  a wire. 
For  this  we  need  to  incorporate  the  heat  gained  by  the  reactant  gas 
stream.  This  forces  us  to  consider  an  energy  balance  for  the  reactant 
gas  stream.  Following  Edwards  et  al.  (1973)»  we  write  the  dimensionless 
mass  and  energy  balance  equations  as 


dc/  _w 

dt  ^ ^A  ■ ^A® 


(C.1  ) 


3t'^  _ 13 V w-  ^ w w^w  T 

3t  - ^ ^ B Da  C^e 

G 0 X 


,W 


(C.2) 


subject  to 


= 0 at  X = 1 

3t'^  w 

- h(T  -T  . ) = 0 at  X = 0 
3x  gin 

In  the  writing  of  these  equations  we  have  neglected  diffusion  of 
active  sites  and  assumed  the  wire  length  is  finite.  The  energy  balance 
equation  for  the  gas  stream  is 
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dT 


Le 


dt  2^^ 


-T  ) + ) 

gin  g g 


(C.3) 


The  superscript  w has  been  used  to  denote  the  variables  and  param- 
eters for  the  reaction  on  a wire  which  are  analogous  to  the  CSTR.  Here 
w . ^ . 

E is  a parameter  inversely  proportional  to  the  thermal  conductivity  of 
the  wire.  We  seek  the  limiting  form  of  these  equations  as  e'^  ->•  0.  The 
state  variables  are  expanded  in  a power  series  of  e'^  as 


t”  = E T':'(e'^)^ 
i=0  ^ 


(C.4) 


Substituting  (C.4)  into  (C.2)  we  obtain 


dc’J  t” 

- 1 - C“  - Da“c“  e ° 
dt  Ao  Ao 


(C.5) 


Le 


dT 

w go 
dt 


ot^(T  . -T  ) + (t'^-T  ) 

2 gin  go  o go 


(C.6) 


8^t'^ 
0 

9x^ 


= 0 


9t”(0)  9t'^(1) 

0 o 


9x 


9x 


(C.7) 


Equation  C.7  yields  T^  is  independent  of  x,  i.e.,  it  is  a constant 
throughout  the  wire.  Its  variation  with  time  is  found  from  the  next 
order  term 


dl” 

-S- 1 " -t")  * B”Da"c“  e “ 

dt  „ 2 go  o Ao 

oX 


9T 


,w 

^ = h(T'^-T  . ) at  X = 0 
9x  o gin 


9T 


,w 


^ = -h(T'^-T  . ) at  X = 1 
9x  0 gin 


(C.8) 
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Integrating 


(C.8)  over  x we  obtain 


dT 
o 

dt 


2h(T  . -t'^)  + (T  -t'^)  4-  e'^Da'^c'^ 


gin 


go 


Ao 


(C.9) 


The  analogy  between 


equations  (C.5, 


C.6,  C.9)  to  equations  (2. 1-2.3)  can 


be  seen  directly. 
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